AN ACCURATE FINITE ELEMENT METHOD FOR ELLIPTIC 
INTERFACE PROBLEMS 



GUNTHER H. PEICHL * AND RACHID TOUZANI t 

Abstract. A finite element method for elliptic problems with discontinuous coefficients is pre- 
sented. The discontinuity is assumed to take place along a closed smooth curve. The proposed 
- - ■ method allows to deal with meshes that are not adapted to the discontinuity line. The (nonconform- 

ing) finite element space is enriched with local basis functions. We prove an optimal convergence 
' rate in the H^-novm. Numerical tests confirm the theoretical results. 

o 

(N . 

1. Introduction. Boundary value problems with discontinuous coefficients con- 
stitute a prototype of various problems in heat transfer and continuum mechanics 
where heterogeneous media are involved. The numerical solution of such problems 
requires much care since their solution does not generally enjoy enough smoothness 
properties required to obtain optimal convergence rates. Although fitted or adapted 
meshes can handle such difficulties, these solution strategies become expensive if the 
\ discontinuity front evolves with time or within an iterative process. Such a (weak) 

y-^ • singularity appears also in the numerical solution of other types of problems like free 

I boundary problems when they are formulated for a fixed mesh or for fictitious domain 

• methods. 

, We address, in this paper, a new finite element approximation of a model elliptic 

transmission problem that allows nonfitted meshes. It is well known that the stan- 
dard finite element approximation of such a problem does not converge with a first 
order rate in the _ff ^-norm in the general case. We propose a method that converges 
optimally provided the interface curve is a sufficiently smooth curve. Our method is 
' based on a local enrichment of the finite element space in the elements intersected by 

If^ , the interface. The local feature is ensured by the use of a hybrid approximation. A 

' Lagrange multiplier enables to recover the conformity of the approximation. The de- 

rived method appears then rather as a local modification of the equations of interface 
' elements than a modification of the linear system of equations. This property ensures 

I that the structure of the matrix of the linear system is not affected by the enrichment. 

Let us mention other authors who addressed this topic in the finite element con- 
text. We point out the so-called XFEM (eXtended Finite Element Methods) devel- 
oped in Belytschko et al. [3j where the finite element space is modified in interface 
►v>( \ elements by using the level set function associated to the interface. Such methods, 

5^ ' that are used also for crack propagation, have in our point of view, the drawback of 

resulting in a variable matrix structure. Moreover, although no theoretical analysis 
is available, numerical experiments show that they are not optimal in terms of accu- 
racy. Other authors like Hansbo et al. [13l [12] , have similar approaches to ours but 
here also the proposed method seems to modify the matrix structure by enriching 
the finite element. In Lamichhane-Wohlmuth [16] and Braess-Dahmen a simi- 
lar Lagrange multiplier approach is used for a mortar finite element formulation of 
a domain decomposition method. Finally, in a work by Li et al [15j . an immersed 
interface technique, inspired from finite difference schemes, is adapted to the finite el- 
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ement context. Note also that the references where Lagrange multiphers are employed 
have for these multipliers as supports the edges defining the interface. In our method, 
the interface supports the added degrees of freedom but the Lagrange multipliers are 
defined on the edges intersected by the interface and thus serve to compensate the 
nonconformity of the finite element space rather than enforcing interface conditions, 
which are being naturally ensured by the variational formulation. 

In the following, we use the space L^(17) equipped with the norm || • ||o,o and the 
Sobolev spaces H^{il) and W'^'P{il) endowed with the norms || • |j„i,o and || • ||m,p,f2 
respectively. We shall also use the semi-norm | • |i^o of H^{U). Moreover, if Sli and 
form a partition of fi, i.e., fl = niUil.2, riinri2 = and if is a function in W™~'^'P{iT) 
with v\n, e VV^^Pin^), then we shall adopt the convention v £ U f^2) and 

denote by ||t'||m,p,siiua2 the broken Sobolev norm 

Ik||m,p,niun2 = ll^^l|m-l,p,0 + lkl|m,p,Oi + 1 1 ^' || m,p,a2 ■ 

Similarly, we denote by || • ||m,f2iun2 and | • |m,aiua2; the broken Sobolev norm and 
semi- norm respectively for the TJ^-space. Finally, we shall denote by C, Ci, C2, . . . 
various generic constants that do not depend on mesh parameters and by \A\ the 
Lebesgue measure of a set A and by A° the interior of a set A. 

Let f2 denote a domain in with smooth boundary F and let 7 stand for a closed 
C^-curve in f2 which separates fl into two disjoint subdomains fl~ such that 
= r2+ U 7 U fi" and dQ+ = 7. For given / e L'^{Q) and a € we consider the 

transmission problem: 



(aV?i) 


= / 


in S1+ 


u 


= 


on F, 




= 


on 7, 



where [v] denotes the jump of a quantity v across the interface 7 and n is the normal 
unit vector to 7 pointing into fl^ . For definiteness we let [v] = 1;^ — v+ with = W| j2± • 
In addition to boundedncss of the diffusion coefficient we assume 

(1.1) 

a(x) > a > 0, for a; e ri, 

i.e. a is uniformly continuous on \ 7, but discontinuous across 7. 

The standard variational formulation of this problem consists in seeking u G 
Hlip.) such that 

I a\7u-\7vdx= [ fvdx V v e H^{n). (1.2) 
Jn Jn 

In view of the ellipticity condition (|l.ip . Problem (jl.2p has a unique solution u in 
HQ{n) but clearly u ^ H'^{Vl). We shall assume throughout this paper the regularity 
properties: 

ll"ll2,o-uo+ <C||/||o.o. (1.3) 
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Note that these assumptions are satisfied in the case where a^fi- and a\fi+ are constants 
(see [MlllH] for instance). 

In the following, we describe a fitted finite element method, defined by adding 
extra unknowns on the interface 7. It turns out that this method leads to an optimal 
convergence rate. Although it is well suited for the model problem it seems to be 
inefficient in more elaborate problems which, for example, involve moving interfaces. 
To circumvent this difficulty, we define a new method where the added degrees of 
freedom have local supports and then yield a nonconforming finite element method. 
We show that the use of a Lagrange multiplier removes this nonconformity and ensures 
an optimal convergence rate. 

2. A fitted finite element method. Assume that the domain is a convex 
polygon and consider a regular triangulation of with closed triangles whose 
edges have lengths < h. We assume that h is small enough so that for each triangle 
T e =3^ only the following cases have to be considered: 

1) Tn7 = 0. 

2) T n 7 is an edge or a vertex of T. 

3) 7 intersects two difi'erent edges of T in two distinct points different from the 



where Pi{T) is the space of affine functions on T. A finite element approximation of 
(jl.2p consists in computing Uh € Vh such that 



It is well known that, since u ^ H'^{Q.), the classical error estimates (see P) do not 
hold any more even though we still have the convergence result. 



A fitted treatment of the interface 7 can however improve this result. Let for this 
purpose 3/'^ denote the set of triangles that intersect the interface 7 corresponding 
to cases 3) and 4) above. 



and consider a continuous piecewise linear interpolation of 7, denoted by 7^, as shown 
in Figure [^m Clearly, 7^ is the line that intersects 7 at two edges of any triangle that 
contains 7. Unless the intersection of 7 with the boundary of a triangle T does not 
coincide with an edge, T is split into two sets and T~ separated by the curve 7. In 
case 3), the straight line 7?; HT splits T into a triangle Ki and a quadrilateral that we 
split into two subtriangles if 2 and i^a, where we choose K2 such that Ki n K2 ~ Ih- 
In case 4), 7^1 fl T splits T into two triangles Ki and K2- In this case we set = 0. 
This construction defines the new fitted finite element mesh of the domain (see 
Figure 12. ip . The splitting T = Ki U K2 U is not unique but the convergence 
analysis does not depend on it. Let us denote by the set of the three subtriangles 





(2.1) 



lim ||u - UhWisi = 0. 
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of T. Below S'h will stand for the set of all edges of elements and S"^ is the set of all 
edges that are intersected by 7 (or 7^), i.e. 

For each T € S/h, S't \s the set of the three edges of T. The fitted mesh is denoted by 



Fig. 2.1. Subdivision of interface triangles. 



and by Sj^ :— 1J{T; T G =^7}- ^et us finally note that the curve 7^, defines a new 
splitting of O into two subdomains SX^ and il^ where is defined analogously to 
51^ with 7 replaced by ^h- 

Next we construct an approximation of the function a on the elements of '■ 
For this purpose, let be extensions of to such that a"*^ € Such 
extensions exist due to the regularity of 7 (see [T]). Define o/i € T4^-'^'°°(f7) by 



in 

in o;^, 



and denote by ah the piecewise linear interpolant of a on . Hence is continuous 
on U Sl^ and coincides with a on the nodes of ■ addition, the function 
is discontinuous across the line 7^ and satisfies the properties, 

G e (2.2) 

I h 1/1 

||ah||o,oo,o < C ||a||o ,00,0 5 

(2.3) 

a/j > a > a.e. in ?L. (2.4) 
We now define the finite element space 

Xh := {v e C"(n); «|a\s- = 0, G Pi{K) V K e .^.^ , V T e ^^J}. 

Note that we have Wh C HQ{fl). A fitted finite element approximation is defined as 
the follows: 

Find e Wh such that 

„ f (2-5) 

ahVui-^ ■ Vvdx — I fvdx V u G Wh- 
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In order to study the convergence of Problem (|2.5p . we consider the auxihary 
problem: 

Find e Hq{Q) such that 

(2.6) 



tth'^Uh ■ '^v dx = / fvdx V V £ H^{Q). 
1 Jn 



We note that both problems (|2.5p as well as (|2.6[) have a unique solution. The 
regularity properties (|1.3|1 imply u+ e C"(r2'*"), € C°(r2~) and that u+ and u~ 
have a common trace on 7. Therefore w is continuous on Q and the piecewise Pi 
interpolant IhU G Wh is well defined. In the following let e H^{Q) stand for the 
extensions of from Sl^ to fi. 

In the sequel, we assume that the fitted family of meshes U ''^h)h satisfies the 
condition 

-<Ch-' (2.7) 
Q 

for some 9 G [0, 1) and for which C is independent of where g denotes the radius 
of the largest ball contained in any triangle in any triangle T G S^]^ ■ 
Lemma 2.1. Let u e H^{n+ U fl^). 
1. We have the local interpolation error 



Ch\u\2,T forT^,%\,^;/ 
C^(\u+\2,K + \u-\2,k) for K e 3r^, T e 



where qk is the radius of the inscribed circle of K. 
2. The global interpolation error is given by 

\u - hu\i,n < Ch^-'' |u|2,n+uo-- (2-9) 

Moreover, if u e W'^'^{Q.+ U Q.^) then 

\u - Ihu\i.Q. <Ch |u|2,oo,n+un- • (2-10) 

Proof. Since the local interpolation error estimate for T & ,'7^ is classic in 
finite element theory (see |5] or [S] for instance), we only need to prove the second 
estimate on triangles where u is only piecewise smooth. Consider an element T € S/'^ 
and any subtriangle K e S^^ . Without loss of generality we assume K C fi^, then 

k = {k^vl-^)\j {Knn-). 

Since K n i}^ C T n fi^ n and jh interpolates the interface 7 we obtain for the 
measure of A' n 

\Knn-\<\Tnn- nn+\<ch^, (2.11) 

with a constant C > which depends on 7 only. In view of IhU — IhU'^ , the standard 
interpolation theory (see [S] or [5]) implies 



\u - Ihu\i^K <\u- u'^\i,K + - IhU'^l 



l.K 



(2.12) 

< \u-U+\i,K+C \U+\2,K. 

QK 



Since u'^ — u holds on if n f7+ we obtain 

\u-u+\i^K = \u- u+\i,Knn- < \u~\i.,Knn- + \u'^\i,Knn- ■ 

Applying Holder's inequality with p — ^ and 9 = 3, the imbedding of H^{K) into 
L^{K) (Note that the imbedding constant can be bounded independently of h) and 
(|2.11|) one can bound \u^\i^Kr\n- (and analogously \u'^\i^Kr\n-) by 

\u^\i,Knn- < \K nn^\^\\\7u-\\o,eMnn- 

< C h\\Vii-\\o,e,K < C h\u-\2,K- 

Hence 

\u - U+\iM < Ch{\il-\2M + \u^kK)- 

Inserting this estimate into (|2.12p leads to 

\u-Ihu\i^K<C — (|m \2,k + \u'^\2,k)- 
Qk 

To prove the global interpolation error bound, we write 

\u-Ihu\\n= ^ \u~Ihu\lT+ ^ ^ \u-Ihu\l^K 

<ch' E Kt + cY. —(\^'\Ik + \u+\Ik) 

<C-i\u-\l,, + \U+\l^) 
g 

h'^ 2 
< C I"l2j2+un-' 



where 







g = mm{gK : K e ,T e ^7}. 



The calculation above indicates how the convergence rate can be improved in case 
u e W^'°°(n+ U n-) observing that \Sl\<Ch holds. □ 

Remark 2.1. It is classic infinite element theory to assume that the meshes are 
regular in the sense that Condition (12. 7p is satisfied for 9 = 0. For the fitted meshes 
S^^ one cannot guarantee that such a condition is satisfied. To relax this constraint, 
we assume here (|2.7p for a 9 [0, 1) thus allowing a larger class of fitted meshes than 
permitted by 9 = 0. 

The following result gives the convergence rate for Problem (|2.5[) . 

Theorem 2.1. Assume that the family of fitted meshes {^if)h satisfies the reg- 
ularity property (j2.7p . Then we have the error estimate 

p ^ jch^-'\\u\\2,in^n- ifueH^{n+un-), 

u — U;, 1 n < s ^ . 2.13 

' [C/ih||2,oo,o+uo- ifu^W^^^{^+^^l-). ^ ' 
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Proof. We have from the triangle inequality 

\u - u^li.n <\u~ + \uh - 



(2.14) 



To bound the first term on the right-hand side of (|2.14p . we proceed as follows: Let 
us subtract (|2.6p from p.2p and choose v — u ~ U}^. We have 



/ (a Vm — ah^Uh) ■ V(m — Uh) dx = 0. 



Then 



a}i\S/{u - {t;i)p dx 



u — Uh) dx. 



= — (a — ah) Vu • V(u — Uh) dx — / (a ~ Oft.) ' V( 

The usual estimate for the interpolation error gives 

||a- a;j||o,oo,o < C'/i(l|a|li.oo.f1+ + ll«lll,oo,f2-) 
< Ch ||a||i^oo,o+uf2-- 

with a constant C which only depends on a reference triangle, (see [8], p. 124). Thus 
we obtain 



{a~ah) Vii- V(ii — w/i) dx 



n\sl 



< C /i|la|li^^^o+uf2- |w|i,n\s;^ \u-Uh\i,n\sl- (2-15) 



Next we consider a triangle T E which we split as 

T = (T n r2+ n 17+) u (T n n ^i) u (r n f7+ n ^i) u (r n f^- n r2+). 

As before, we obtain 



Tnsi+no+ 



(a— ah) Vu-V(m— u/i) dx 



< C/i ||a|li^oo^j^+un- l"li,Tnn+nn+ l"~'"''li,Tno+nn+- 



Arguing as in the proof of Lemma[2TTJ the generalized Holder inequality together with 
(|2.1ip yields the estimate 



(a — a/i) Vu • V(u — Uh) dx 



Tnn+nn^ 

= / (a+ — a^)VM+ • V(w+ — w/i) c?x 

< C llallo.oo.o |Tn r!+ n W^u+W^ ,^^^^^^^^- ||V(w+ - Uh)\\o^Tnn+nn- 

<Ch ||a||o,oo,o |w^|2,T ||V(m+ - M/i)||o,T- 

Analogous estimates hold with + and — interchanged. Collecting the four contribu- 
tions to the triangle T one obtains 



(a — ah) Vu ■ V (u ~ Uh) dx 

< Cft.(||a||o,oo,o + ||a||i,oo,f2+ua-) 

X (|u+|2,T ||V(u+ ~Uh) ||o.Tnn+ + \u^\2,T ||V(m" -Uh)\\o,Tnn-)- 

7 



Combining this estimate with (|2.15p leads to 

J ah |V(u - uh)f dx < C/i|la|li_oo,n+ua- |M|i,n\s,-^k - "/i|i,o\s^ 
+ Ch (||a||o,oo,fi + ||«||i,oo,o+ufi-) 

X X! (|"^|2,t||V(u+ - u/i)||o,Tno+ + |w"|2,t||V(w" - u/i)||o,Tnn-) 

<Ch ||a|li,oo,n+un- l"li,o\s,; I" - ^^ft|i,n\s,"[ 

+ Ch{\\a\\a^oosi. + \\a\\i,oo,n+yjn-)(\u^\2,sl + l^^^.s,-:) - Uh)\\a^sl 
< C'/i (||a||o,oo,J2 + ||a||i^oo,fi+uJ2-) |u|2,o+uf2-||V(M- u/i)|lo,n, 

which by ()2.4p impUes 

I" - uh\i,n < C/i (||a||o,oo,f2 + ||a||i,oo,o+uJ2-) Iwb^n+ua-- (2-16) 
To bound the norm \uh — u^\i,n, we have from problems (|2.6p and (|2.5p . 

Standard finite element approximation theory combined with (|2.2p - (|2.3p gives 

\uh-ul\i^n<C inf \uh~ v\i^q,, (2-17) 

which together with (|2.16p implies 

\uh - |i,o < C \uh - //iu|i,n 

< C \uh - + C |w - //iu|i,n 

< C/i |u|2,o+uo- +C|u-/?,M|i,n- 

The interpolation error is bounded using (|2.9p or (|2.10p . □ 

3. A hybrid approximation. The method presented in the previous section 
has proven its efficiency as numerical tests will show in the last section. In more 
elaborate problems like time dependent or nonlinear problems where the interface 
7 is a moving front, the subtriangulation moves within iterations and then the 
matrix structure has to be frequently modified. To remedy to this difficulty, we resort 
to a hybridization of the added unknowns. More specifically, the added discrete space 
Xh is replaced by a nonconforming approximation space. In addition, a Lagrange 
multiplier is used to compensate this inconsistency. The hybridization enables to 
locally eliminate the added unknowns in each triang le T G In the sequel we 

fix an orientation for the interface 7. This induces an orientation of the normals to 
the edges e S S'^ by following the interface in the positive direction. The jump of a 
function v across an edge e G S'^ can then be defined as 

[u]e(a;) '■— lini v{x + sn{x)) — lini v{x + sn{x)) = v'^ [x) — {x), x € e, 

s^O,s>0 s^O,s<0 

where n is the unit normal to e. 
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To develop this method, we start by defining an ad-hoc formulation for the solu- 
tion Ufi of (|2.6p . Let us define the spaces 

Yh {v e L\n)- v\n\si = «|t e H\T) V T e 
[v] = on e, V e e <4 \ 4^}, 

where -ffoo^ (^) ^^^^^ space of the trace space 

ij|,(e) := t; e i?i(r), e G fSr, u = on d V d G ,#t, 7^ e}. 

^ i 
We remark that the jumps [v\ for w G Zh can be interpreted in iJgQ(e) for e G f?^. 

This is due to the fact that v G H'^{T) for all T e that for every e G the 

jump of V lies in H'2{e) and vanishes at the endpoints of e as well as on at least two 

adjacent edges. This motivates the choice of Qh- 

The elements of Qh will be referred to by = (/^e)ee^J- We endow Zh with the 

broken norm 

ii"iiz. = (E \<tY''- 



On Qh we use the norm 



^C.^'T -"00 '.'^^ \ 



^eW ds 



sup 



2\ 2 



Above, the integrals over edges e are to be interpreted as duality pairings between 

— 11 ^ 
iJgQ^ (e) and i?Q^Q(e). We mention that the broken norm in Zh reflects the fact that 

Zh is not a subspace of Hq{VI). 

Next we define the variational problem. 

Find (u^. Ah) G Zh x Qh such that: 
^ / a;i Vuf • Vw dx - ^ \h[v]ds^ / /wdx V w G (3.1) 

^ [fi[uj^]ds = y^eQh- (3.2) 



The saddle point problem (|3.ip - p.2p indicates that the continuity of m/i across the 
edges of S'^ is enforced by a Lagrange multiplier technique. 

Theorem 3.1. Problem p.ip - p.2p /las a unique solution {u^ , \h) G x Q/j. 
Moreover, we have = Uh and the following estimate holds 

\\uUz, + \\\h\\Q^<C\\fh,n, (3.3) 



with a constant C which is independent of h. 

Proof. Problem (|3.ip - (|3.2p can be put in the standard variational form 



where 



fv dx. 



The bilinear form is clearly continuous and coercive on the space Zh x Zh. The 
bilinear form ^ is also continuous on Zh x Qh- 

Next we verify that !3§ satisfies the inf-sup condition, i.e. there exists 5 > Q such 
that for every X (z Qh there exists £ Zh such that 



I.e. 



holds. 



J2 / t^e[v^,]ds >S\\Vf,\\gJ\fl\\Q^ 



(3.4) 



Given ^ = (Me)ee<f^ ^ Qh and an edge e G S'j^ choose a triangle T S which 
has e as one of its edges. Define vt G H^{T) as the solution of 




in T, 
on e, 

on ar \ e, 



(3.5) 



which is equivalent to 



j \7v-\7ipdx = J iJe'fds for e iJg(T) 



where 

Hl{T) ^{ipe H\T); V? = on dT\e}. 
By Green's theorem we obtain 

dvT 



\\t^e\\-l/2,e 



dn 



-1/2,, 
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which impHes 

||A^e||^l/2,e < J \Vvt\'^ dx ^ J^iieVrds. 
Let XT denote the characteristic function of T and define 

Y XtVt- 



Since there are as many edges in as triangles in then [v^] — vt holds for every 
edge e S (f^. Hence we obtain 

Furthermore, 

II^mIII,- Y W^^tWIt^ Y W^^tWIt 

holds. This implies 

ii/^iilJi^/^iii ^ ( E w^vtWItY -mv^,f^r■ 
Adjusting the sign of Vf_i this is equivalent to (|3.4p with 6 = 1. The estimate p.Sp is 
a direct consequence of p.4p . 

Now, it is clear from p.2p that 

[wf ] =0 on e, V e e . 

This implies that uf^ S iJQ(Sl). Choosing a test function v E -ffo(f^) in (13. ip . we find 
that is a solution to Problem (|2.6p . and then = w^i. The interpretation of 
is simply obtained by the Green's formula. □ 

We are now able to present a numerical method to solve the interface problem. 
This one is simply derived as a finite element method to solve the saddle point prob- 
lem p.ip - p.2p . We consider for this end a piecewise constant approximation of the 
Lagrange multiplier. Let us define the finite dimensional spaces, 

Zh ■.= Vh + Yh, 

Yh := {v e L\n); v\n\sZ = 0- ^'1^ ^ Pi{K) y K e .^^ , ^ T € ^J, 

[v] = on e, V e e \ ^^^}, 
Qh := {/i e Yl ^^(e); f^\e = const. V e e 4^}. 

The hybrid finite element approximation is given by the following problem: 
Find (u^ , Xh) G Zh x Qh such that: 

"Y^ / o.h'^u^ ■ '^v dx ~ "Y^ / ^h[v]ds= / fvdx y v G Zh, (3.6) 
^ l'^i[u^]ds = ytieQh- (3.7) 

11 



Let us give some additional remarks before proving convergence properties of this 
method. 

1. The matrix formulation of the method has the following form 




ACQ" 

B^\\v\ = \c\, (3.8) 

where the vector u contains the values of at nodes of the mesh i.e. components 
of in the Lagrange basis of V/i, v contains the components of in the basis of 
Y/i, and A has as components the values of A/i on the edges of S'^. There is clearly 
no simple method to eliminate off diagonal blocks in the system p.8p in order to 
decouple the variables. More specifically, our aim is to eliminate the unknowns v. 

2. The method must be viewed in the context of an iterative process like the 
Uzawa method, where the Lagrange multiplier A^ is decoupled from the primal vari- 
able . In such situations, each iteration step consists in solving an elliptic problem 
with a given A^. Let us recall that, due to the local feature of the basis functions of 
nodes on edges of S'^ , the unknowns associated to these nodes can be eliminated at 
the element level. This is a basic issue in our method. 

3. We point out that equation (|3.7[) entails 

[uf ] =0 on e, V e e 5^'^. (3.9) 

This follows from the fact that is an afhne function on each edge of . This 
implies that actually G Wh. Choosing v £ Wh in (|3.ip we find 



/ aiiVuf^ ■ Vv dx — / fv dx. 
Jn Jn 



This yields uf^ = u^. 



4. Convergence analysis. This section is devoted to the proof of existence, 
uniqueness and stability of the solution of l|3.6p - p.7p as well as its convergence to 
Problem (IXB^d^. 

For this result we need a localized quasi-uniformity of the mesh. More precisely, 
we assume that 

\e\>Ch V e e . (4.1) 

In addition, we make the following assumption: 

The distance of the intersection point of 7 with any edge e € S"^ 

to the endpoints of e can be bounded from below by Sh, where S is (4-2) 

independent of h. 

Although this assumption appears to be quite restrictive, numerical tests have shown 
that it can be actually ignored in applications. 

Theorem 4.1. Assume that the family of meshes {,'7h)h satisfies Property (j4.ip . 
Then Problem (|3.6p - (l3.7p has a unique solution. Moreover, we have the bound 

\\u?Az, + \\\h\\Q,<C\\fh,n. (4.3) 
12 



where the constant C is independent of h. 

Proof. It is clearly sufficient to prove the inf-sup condition (see for instance 
Brezzi-Fortin [7]): 

sup „^ ;; \, I, >/3>0 V/i,.eQ/.. (4.4) 

In the following, for each triangle T e , we shall denote by (resp. e^) the edge 
where 7 enters T {resp. leaves T), and by ex the remaining edge of T (see Figure 
(|4.ip ). Recall that we fixed an orientation for 7. 



ih 



T+ 



Fig. 4.1. Definition of e^, e^,, and ex- 



Let /i/i e Q/i, and let u e Z/j be the function given by Problem p.Sp . 
We define a function G Zh by 

«^|T = V T e \ 

Vhds^ I vds \/ eec^T, y T e 



(4.5) 



The gradient of can be expressed in T € .^f^ by 



where 95^+ {resp. ip^- ) is the basis function of Zh associated to the added node on 
{resp. eTp). Then by using (14. ip and the Cauchy-Schwarz inequality, we get for each 



|Vv,,||o,T = Cih- 



< ds 



V(/3g- ||o,T + C2 h' 



' ds 



< C3 h'-- {\\v\\o,e- II V^e- llo.T + |k||o,e+ II V^^.J ||o,t). 



(4.6) 



The trace inequality (see [1], eq. (2.5)) and the Poincare inequality owing to w = 
on ex, yield for T e ^J, 

lkllo,e± < Ci {h--2\\v\\o,T + h^\\yv\\o,T) <C5/i^||Vi;|lo,T. (4.7) 
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On the other hand, Assumption (j4.2p imphes the uniform boundedness of || Viy9g± ||o,t- 
From (|4.6p and (|4.7p we obtain then 

||Vwft||o,T < Cq ||Vw||o,t- 
Using the inf-sup condition p.4p and (|4.5p . we finaUy obtain 

II^'^HqJI^'^IIz. - C'ellAi/JIgJK'b, 

< C7 ^ / [w] ds 



C7 ^ / K] ds. 



FinaUy, obtaining the estimate (|4.3p is a classical task that we skip here. □ 
We now prove the main convergence result. 

Theorem 4.2. Assume hypotheses p.7p and (|4.ip are satisfied, then there exists 
a constant C , independent of h, such that 



Proof. From classical theory of saddle point problems (see [T^, p. 114), we obtain 
from Theorem 14. 11 



- ^hWz, + - Xh\\Q, < C ( m|^ - v\\g^ + ||A, - mHqJ. (4.8) 

Furthermore, using (Braess [4], Theorem 4.8), Property p.9p implies that Estimate 
(|4.8p can be improved, for the error on uf^ by 

\\uf^ -uj^\\z,<C inf ll^lf (4.9) 

To bound the right-hand side, we choose v — luu, where Ih is the previously defined 
Lagrange intcrpolant in Z/,. Since = Uh (see Theorem 13. ip . then by using (|2.9p 

and (mm), 

< \\u-Ihu\\2^ + \\u-Uh\\z^ 

< Ci h^^^ |'u|2,n+un- + C2 /i |u|2,n+un- • 
If u G W'^'°°{Sl+ U rj-), then Estimate (piU)) yields 

< C/i||u||2,oo,n+un-- (4.10) 

□ 

Remark 4.1. As it was previously mentioned, we know that if Problem p.ip - 
()3.2p has a unique solution {u^ , \h) then uf^ — where is the solution of Problem 



and therefore the error estimate (I2.13P holds. Consequently, Theorem \4.S\ can be 
simply proven by obtaining a nonuniform inf-sup condition (i.e. (14. 4p with [3 — (3(h)). 
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This can be achieved without assuming (j4.2p . In this case, no error estimate is to he 
expected for the Lagrange multiplier. 

Finally, since the Lagrange multiplier can be interpreted in terms of Uh (see 
Theorem 13. it is interesting to see how good is its approximation Xh- Let, for this. 
Eh denote the set 



Eh-.^ n 



e. 



Theorem 4.3. Under the same hypotheses as in Theorem \4-^ we have the 
following error bounds 



'C{h^-' + hi) {\u\^^n+un- + l|A||o.B.) ^fXh e L^iE^), 
Ch^-' (|u|2,n+uo- + m^.E,) ^f^h e Hi{Eh). 



Proof. We use the abstract error bound (|4.8p . Let, for e E S'^ , 

Ae := 7^ [Xhds. 

|e| Je 

Using Lemma 7 in Girault-Glowinski we obtain the bound 

Uh-Xe\\ 1 <C;i^||A,,||o.e iiXhEL^ie), 

H„o (e) 

and 

||A^-Ae|L -1 <Ch\\Xh\\i,, ifX^eHHe), 

Combining these bounds with (|4.8p . (|4.9p and (|4.10p achieves the proof. □ 

5. A numerical test. To test the efficiency and accuracy of our method, we 
present in this section a numerical test. We consider an exact radial solution and test 
convergence rates in various norms. 

Let il denote the square ft ~ (—1,1)^ and let the function a be given by 



a{x) 



a if |a;| < 
P if \x\ > Ri, 



where a,f3 > 0. We test the exact solution 

U{x) : 



^iRl-\xn + ^iRl-Rl) if N<i?i, 
j^iRl-\x\') if N>i?i. 



We choose i?i = 0.5 and i?2 = V^- The function / and Dirichlet boundary conditions 
are determined according to this choice. Note that unlike the presented model prob- 
lem, we deal here with non homogeneous boundary conditions but this cannot affect 
the obtained results. 



15 



The finite element mesh is made of 2N'^ equal triangles. According to the defini- 
tion of a, the interface 7 is given by the circle of center and radius The error is 
measured in the following discrete norms: 



|e||o,/i 



1 

M 



M 



i=l 



|e||o,oo := max \u{xi) ~ Uh{xi)\, 

l<i<M 



/ \Ih{\/u){x)-'^Uh 



where Xi are the mesh nodes, M is the total number of nodes, and Ih is the piecewise 
linear interpolant. We denote in the sequel by p the ratio a/f3. Table 1 presents 
convergence rates for the standard Pi finite element method using the unfitted mesh 
?T|) with the choice p = 1/10. 





l|e||o,h 


Rate 


l|e||o,oc 


Rate 


Mi.h 


Rate 


10 


1.40 X 10-^ 




2.02 X 10"^ 




6.28 X 10-^ 




20 


6.78 X 10-^ 


1.05 


1.09 X 10"^ 


0.89 


5.23 X 10"^ 


0.26 


40 


3.61 X 10"^ 


0.91 


5.81 X 10"^ 


0.91 


3.68 X 10-^ 


0.51 


80 


1.83 X 10-^ 


0.98 


3.06 X 10"^ 


0.92 


2.56 X 10"^ 


0.52 


160 


9.44 X 10-^ 


0.95 


1.55 X 10"'^ 


0.98 


1.82 X 10-^ 


0.49 



Table 1. Convergence rates for a standard (unfitted) finite element method. 

As expected, numerical experiments show poor convergence behavior. Let us con- 
sider now the results obtained by the present method, i.e. p.ip - (|3.2p or equivalently 
(|2.5p . We obtain for p — 1/10 and p = 1/100 the convergence rates illustrated in 
Tables 1 and 2 respectively. 





\\e\\o,h 


Rate 


e 0,0c 


Rate 




Rate 


10 


3.45 X 10-^ 




4.25 X 10-=* 




1.75 X 10-^ 




20 


8.18 X 10-'' 


2.1 


1.72 X 10"^ 


1.3 


6.87 X lO--' 


1.3 


40 


1.70 X 10-* 


2.3 


5.22 X 10-'' 


1.7 


2.81 X lO"-^ 


1.3 


80 


3.94 X 10-^ 


2.1 


1.64 X 10"* 


1.7 


1.02 X 10-=* 


1.5 


160 


8.57 X lO-t- 


2.2 


4.89 X 10"^ 


1.7 


3.59 X 10-"* 


1.5 



Table 2. Convergence rates for the hybrid finite element method with p = 1/10. 





M\o,h 


Rate 


e 0,0c 


Rate 


l|e||iJi 


Rate 


10 


3.26 X 10-^ 




4.07 X 10-^ 




1.69 X 10-^ 




20 


7.91 X lO-'' 


2.0 


1.74 X 10"=* 


1.2 


6.65 X 10"^ 


1.3 


40 


1.72 X lO-'' 


2.2 


5.47 X 10-* 


1.7 


2.72 X 10-^ 


1.3 


80 


4.01 X 10-^ 


2.1 


1.74 X 10-'' 


1.6 


9.88 X 10-" 


1.5 


160 


8.82 X lO"*^ 


2.2 


5.22 X 10-^ 


1.7 


3.50 X 10-* 


1.5 



Table 3. Convergence rates for the hybrid finite element method with p = 1/100. 

Tables 2 and 3 show convergence rates that are even better than the theoretical 
results. This is probably due to the choice of a discrete norm but may also be due to a 
superconvergence phenomenon. Rates for the L^-norm give also good behavior. For 
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the //""-convergence rate, we can note that we dot retrieve the second order obtained 
for a continuous coefficient problem. However, these rates are better (1.5 rather than 
1) than the ones obtained for a standard finite element method and, moreover, the 
error values are significantly lower in our case. It is in addition remarkable that the 
error values depend very weekly on p but the convergence rates are independent of 
this value. 

6. Concluding remarks. We have presented an optimal rate finite element 
method to solve interface problems with unfitted meshes. The main advantage of 
the method is that the added unknowns that deal with the interface; singularity do 
not modify the matrix structure. This feature enables using the method in more 
complex situations like in problems with moving interfaces. The price to pay for this 
is the use of a Lagrange multiplier that adds an unknown on each edge that cuts 
the interface. This drawback can be easily removed by using an iterative method 
such as the classical Uzawa method or more elaborate methods like the Conjugate 
Gradient. The good properties of the obtained saddle point problem enable choosing 
among a wide variety of dedicated methods. This topic will be addressed in a future 
work. Let us also mention that the present finite element method does not specifically 
address problems with large jumps in the coefficients. These ones are in addition ill 
conditioned and this drawback is not removed by this technique. 
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